A semiquantal approach to finite systems of interacting particles 
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A novel approach is suggested for the statistical description of quantum systems of interacting 
particles. The key point of this approach is that a typical eigenstate in the energy representation 
(shape of eigenstates, SE) has a well defined classical analog which can be easily obtained from the 
classical equations of motion. Therefore, the occupation numbers for single-particle states can be 
represented as a convolution of the classical SE with the quantum occupation number operator for 
non-interacting particles. The latter takes into account the wavefunctions symmetry and depends on 
the unperturbed energy spectrum only. As a result, the distribution of occupation numbers Us can 
be numerically found for a very large number of interacting particles. Using the model of interacting 
spins we demonstrate that this approach gives a correct description of Us even in a deep quantum 
region with few single-particle orbitals. 
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PACS numbers: 05.45.-a 

In many physical systems such as complex atoms, 
heavy nuclei and interacting spins, highly excited eigen- 
states in the unpertured many-particles basis can be 
treated as a chaotic superposition of a very large number 
of components, for a relatively strong interaction among 
particles (see, e.g. This fact has been used in 

in order to develop a statistical description of closed sys- 
tems with a finite number of Fermi-particles. In partic- 
ular, it was analytically shown that for a strong enough 
interaction, a smooth dependence of occupation numbers 
on the energy occurs, which is related to global properties 
of chaotic eigenstates. 

As is known, the direct numerical computation of ex- 
cited eigenstates is strongly restricted due to a very large 
number of many-particles states, which grows extremely 
fast with the number of particles. However, the mean 
values of occupation numbers turn out to depend on the 
average shape of chaotic eigenstates in the unperturbed 
basis, not on exact, specific values of their components 

i- 

In this Letter we develop a novel approach to quan- 
tum systems with chaotic behaviour in the classical limit. 
This approach takes into account both the chaotic prop- 
erties of the classical system, and the specific features of 
the unperturbed single-particle spectrum. As a result, 
one can avoid diagonalization of Hamiltonian matrices of 
a huge size which may be practically unfeasible. This 
kind of approach can be applied to generic Hamiltonian 
systems with two-body interaction of the type: 



H = Ho + V; 



N 

Ho = 



V = 



N 



N 
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V^,r (1) 



Here Hq describes N non-interacting particles with 
as single-particle Hamiltonians, and V stands for a long- 
range two-body interaction between the particles. In 
what follows, we assume that the single-particle spectrum 
is determined by a finite number L of single-particles en- 
ergies eii, h = 1, L; however, the approach is valid for 
more generic systems with an infinite spectrum. 

The unperturbed Hamiltonian Hq determines the 



many-particle states \K) 



|0) (with al 



as creation-annihilation operators), that form the basis 
in which the exact eigenstates of H are represented. As 
usual, we assume that the basis is ordered according to 
increasing unperturbed energy values Eq = X]ft=i ^h- 

The distribution of occupation numbers (DON) for 
single-particle states is defined by the relation. 
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K=l 



where fih 



a\ah is the occupation number operator giv- 



ing the occupation numbers n'^ = {K\fih\K) equal to 
or 1 for Fermi-particles, and to 0,l,2,...,iV for Bose- 
particles. These numbers indicate how many parti- 
cles in a many-particle basis state \K), occupy a partic- 
ular single-particle state \h). Correspondingly, the occu- 
pation numbers UE^h) give the probability that one of 
N particles in a many-particle exact state with the to- 
tal energy E, occupies a particular single-particle state 
\h). The total number M of many-body states equals 
M = L\/{N\{L - N)\) for Fermi and M = {N + L - 
l)!/(iV!(L- 1)!) for Bose-particles. 

One should note that while in the above expression 
for the DON, the eigenfunctions ip^ depend on the to- 
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tal Hamiltonian H, the term nj^ depends on the unper- 
turbed spectrum only. This fact is crucial for our semi- 
quantal approach. Due to the chaotic structure of exact 
eigenstates, one can make an average of the DON over 
a small energy window AE around the fixed value E. 
This averaging procedure is similar to that used in the 
conventional statistical mechanics developed for systems 
with a finite number of particles in contact with a heat 
bath, or for isolated systems of an infinite number of 
non-interacting particles. 

Expression (|^) for the mean values nE{h) can be con- 
siderably simplified by introducing the so-called shape 
of eigenfunctions, SE, (envelope of eigenstates in energy 
representation) . The form of the SE has been studied in 
details both in the model with random two-body inter- 
action pI, and in dynamical models of interacting par- 
ticles ||6|,[7|. The introduction of the average quantity 
SE (thus neglecting correlations between different com- 
ponents V'k) represents the key point of our approach. 

We assume that the unperturbed many-body energy 
spectrum has an intrinsic degeneracy. This situation is 
typical for spin systems, and is more complicated in com- 
parison with those studied before 1^. Below we show 
how this difficulty can be overcome. Let us redefine the 
state \K) by means of a indexes pair |j, ^), where j labels 
the "unperturbed energy" E^ of the many-body state, 
while I — l,...,Nj labels its degeneracy Nj. If there 
are A'o different "unperturbed" energies, one can write, 
T,K=i = Eil'i' therefore, one has. 

No Nj 

nE{h)=Y^J^{j,l\hf,\j,l)\^l;f^f. (3) 

j=i 1=1 

According to standard definition the SE is given 
by: 

N, 

WE{E^)-J2\'J'f.i\'- (4) 
1=1 

By substituting - = We{E^)/Nj with 

< ... >i as an average over I, we obtain an approximate 
expression for the DON in terms of the SE, 

No , 

nE{h) =J2 — Y,{jJ\h^\j,l)WE{E^). (5) 

j=i ^ 1=1 

Needless to say, if an unperturbed spectrum has no de- 
generacy, Eq.(^) can be written in a similar way by taking 
an average over a small window of energy around Ej. 

As one can see, expression (^) depends on two terms 
of different nature. The first one, {j , l\nh\j , I) , refers 
to the unperturbed many-particle spectrum and reflects 
the specific properties of a single-particle spectrum, as 
well as quantum features related to Fermi-Dirac or 
Bose-Einstein statistics. In contrast, the second term, 



We{Ej), refers to global properties of eigenstates and 
describes interaction effects. Therefore, the basic idea 
of our "semiquantal" approach is to substitute the lat- 
ter term (SE) by its classical analog which can be easily 
found from classical equations of motion. 

Classical analogs of the SE have been studied in differ- 
ent models, see, for example, ||J^). In practice, one has 
to derive the distribution WsiEo) = P{Hq = Eq\H = E) 
for the probability to find the unperturbed energy Eq for 
i?o, given the conserved total energy E. This can be ob- 
tained by generating many different initial classical con- 
figurations on the energy surface H ~ E ov sampling the 
i/o(i) values generated by one single trajectory onto the 
energy surface and computing the correspondent distri- 
bution of Hq = i?o §. The two procedures have been 
found to give the same results in the chaotic region [Q. 
In order to facilitate the comparison with the quantum 
SE, in our numerical simulations the bin size of Eq equals 
the energy distance between neighbour values of E'^. In 
the same way one can define the classical distribution of 
occupation numbers, nE{h) = P{h\^ = h\H = E) (see 
also [Ql). For the quantum-classical comparison, the bin 
size of h is taken to be equal to the spacing between close 
single-particle energy levels en- 

Let us stress that in this semiquantal approach (SA) it 
is possible to study specific systems of, let us say, 1000 in- 
teracting particles occupying 10—20 single-particle levels. 
Surely, one expect this approach to be valid for highly ex- 
cited chaotic states in a deep semiclassical region. How- 
ever, by direct numerical simulations we have found that 
the SA gives correct results even for energy values close 
to the ground state. 

The model that we have studied, consists of N 3-D 
interacting spins placed in a magnetic field B directed 
along the z-axis. In order to have a proper many-body 
operator one should require a coupling between all spins 
(not only between neighbors). The Hamiltonian thus 
reads, 

N N~l N 

(6) 

Using the well known relations ~ Sj ± iS^, one can 
write, 

N 

1=1 ] 

N-l N 

^{J. -Jy)E E [S^S+ + Srsr]+ (7) 

i=l j=i+l 
N-1 N 

i=l j=i+l 

In the mean field approximation, we put B+J^ E "^1 — 

j 

const = 1. The interaction can be further simplified by 
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the particular choice Jx = —Jy = J ■ Thus, our Hamilto- Here = aj^ah, with aj^ and ah the creation-annihilation 



nian H = Hq + V has the following form, 



N 



H 



i=l 



N-1 N 
i—1 



(8) 



For simplicity, the N classical constants of motion l^^l 
have been set to 1, therefore, the only free classical pa- 
rameters are the total conserved energy E and the in- 
teraction J. The classical model has been studied in |^] 
and it was numerically found to be chaotic and exponen- 
tially unstable in a wide energy range. More precisely, 
in order to have strong chaos, one needs the interaction 
J between particles to be strong enough. A convenient 
choice is to take the interaction strength J = 1/N. One 
should stress that this situation is the most difficult for 
theoretical studies, see discussion in j?]. 

Quantization follows the standard rules, which gives, 
Sf = h'^m{m + 1) and S'f = hs with — m < s < m. The 
action of creation and annihilation operators is defined 
by: 



.) = h^ym{m+l)- s,{s,±l)\ 



±1,...) 



where |si, sat) are the non-symmetrized states (first 
quantization states). 

There are L ~ 2m + 1 single-particle energy levels 
e/i = —hh with h ~ — to, ...,m. Therefore, the un- 
perturbed many-particle energy spectrum consists of a 
number of degenerate levels with the spacing equal to h. 



Note, that both the ground state Eg 



-mNh and the 



upper level Eu — mNh are non-degenerate. The classi- 
cal limit is recovered when spins are allowed to have any 
possible orientation, that is to — > oo, and ?i — > 0. 

The choice we have done {Jx = —Jy = J) allows to 
reduce the dimension of the Hilbert space by, approx- 
imately, one half. This happens because the operator 
V in (^) connects only those unperturbed many-body 
states that are separated by the spacing 2h. In what 
follows, we consider the subset of the many-body states 
containing the ground state. From these states we con- 
struct the completely symmetrized states |si, .., sat)"^ — 
\n-rm ■■■,nm), where the r.h.s. refers to their second 
quantization representation. Note, that for the sym- 
metrized states the distribution of occupation numbers 
is expected to be described, for a large number of parti- 
cles, by the Bose-Einstein statistics. 

Hamiltonian (|^), in the second quantization, can be 
written as : 



ehnh 



(9) 



with 



h——m 



m-2 

ahah+ X £,hal_^^ah- 

h——m 



operators satisfying the standard relation [a/i,aj,] — Shk- 
As to the coefficients ijh, S,h, they have quite complicated 
expressions; however, they can be easily computed nu- 
merically. 

The procedure we have used in our numerical simu- 
lations consists of the following steps: a) compute the 
classical values nsih), and WEiEo) as described above; 
b) compute the quantum values nsih) and WEiEg) by 
diagonalization of the total Hamiltonian ; c) compute 
nE{h) and We(,Eo) by using the semiquantal approxima- 
tion according to expression (H). 
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(c) 
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FIG. 1. Occupation number distribution obtained in three 
different ways. Crosses correspond to the quantum data, open 
circles stand for the SA approach, and hystograms refer to the 
classical quantities. The data are given for 3 interacting par- 
ticles: (a) E = -2.75, AE = 0.15, m = 5, (h) E = -2.75, 
AE = 0.15, m = 13; (c) E = -0.91 AE = 0.09, m = 5, and 
(d) E = -0.91 AE = 0.09, m = 11. 



The results for the DON are summarized in Figj^. 
For the convenience of comparison, symbols refer to 
the quantum and SA results, while the classical data 
are presented as hystograms,. Note, that each particle 
can have energy within the interval [—1,1]. In order 
to make the quantum-classical comparison as close as 
possible, we took the classical bin size equal to h. All 
distributions have been normalized in such a way that 
J2T=-m.''^E{h) = N. In cases (a) and (b) we choose the 
same energy, close to the ground state and two different 
TO values. In (c) and (d) instead, we choose a higher en- 
ergy value. The latter case( (c) and (d)) should describe 
a situation more classical than the former (higher energy 
should correspond to higher temperature). One can see 
that, while classical and quantum data only disagree in 
the deep quantum region (a) (energy close to the ground 
state and small m (big h)), there is a very nice corre- 
spondence between quantum and SA data in all cases 
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(including region (a)). This nice correspondence is far 
from trivial since SA does not take into account quan- 
tum correlations inside exact eigenstates, together with 
possible correlations between the two terms, {j,l\nh\jil) 
and \tpfi\'^ in expression both expected to be strong 
in the deep quantum region. 

It is very interesting to explore the occurence of the 
Bose-Einstein (BE) distribution in our model. A similar 
problem has been studied in detail for the model of two- 
body random interaction where the conditions for the 
appearence of Fermi-Dirac distribution have been found 
for few interacting Fermi-particles. 

By assuming, a priori, the validity of the BE- 
distribution n^^{i) = [e'^^^'"'^) — l]^-"^ in our closed sys- 
tem, one can find the "temperature" and the "chem- 
ical potential" fi via the standard relations, 

m m 

E n|^(z)=iV; E n^^{{)e,^E'. (10) 

2— — m i— — m 

Here E' is the numerically computed energy obtained 
from the single-particle quantum distribution (see details 
in [^), and N is the number of particles. Notice that, due 
to interaction, E' / E. 



in g). Thus in our model there is no need to increase 
further the number of particles, since even for TV = 10 
a remarkable agreement with BE distribution has been 
found. 

We would like to stress that an exact quantum treat- 
ment of the last example calls for a diagonalization of 
a huge matrix of size 8008 x 8008, while, with the semi- 
quantal approach, all computations required few minutes 
on a standard portable PC. 

In conclusion, we suggest an effective semiquantal ap- 
proach to closed systems of interacting particles, based 
on the chaotic structure of eigenstates. In this approach, 
the computation of the distribution of occupation num- 
bers can be easily performed by making use of the clas- 
sical analog of the shape of eigenstates in the unper- 
turbed many-particle basis. We demonstrate the effec- 
tiveness of this approach using the model of 3-D spins 
with anisotropic Ising interaction. The data show that 
semiquantal computations give results which are very 
close to the exact ones. 
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FIG. 2. Distribution of occupation numbers for N = 10 
particles. Open circles indicates the BE-distribution while 
crosses stand for the SA data: (a) E = —8, AE — 0.2, m — 3, 
(b) E = -6.5, AE = 0.2, m = 13. 

The comparison between BE and SA distributions are 
shown in Fig. 2. As one can see, even for relatively small 
N — 10, the distribution nE{h) can be closely approxi- 
mated by BE-distribution. This confirms the expectation 
that a strong enough interaction between particles can 
play the role of an internal heat bath ||]. Therefore, the 
standard quantum distributions can be used, with a cor- 
responding renormalization of the energy E' (see details 
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